library(knitr)
library(raster) # Raster Data handling## Loading required package: splibrary(tidyverse) # Data Manipulation## ── Attaching packages ────────────────────────────────────────────────────────────── tidyverse 1.3.0 ──## ✓ ggplot2 3.3.2 ✓ purrr 0.3.4
## ✓ tibble 3.0.3 ✓ dplyr 1.0.2
## ✓ tidyr 1.1.1 ✓ stringr 1.4.0
## ✓ readr 1.3.1 ✓ forcats 0.5.0## ── Conflicts ───────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
## x tidyr::extract() masks raster::extract()
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()
## x dplyr::select() masks raster::select()library(getlandsat) # keyless Landsat data (2013-2017)
library(sf) # Vector data processing## Linking to GEOS 3.8.1, GDAL 3.1.1, PROJ 6.3.1library(mapview) # Rapid Interactive visualization## GDAL version >= 3.1.0 | setting mapviewOptions(fgb = TRUE)##Question 1
palo_and_area = read_csv('~/github/geog176A-daily-exercises/uscities.csv') %>%
filter(city == 'Palo') %>%
st_as_sf(coords = c('lng','lat'), crs = 4326) %>%
st_transform(5070) %>%
st_buffer(5000) %>%
st_bbox() %>%
st_as_sfc()## Parsed with column specification:
## cols(
## city = col_character(),
## city_ascii = col_character(),
## state_id = col_character(),
## state_name = col_character(),
## county_fips = col_double(),
## county_name = col_character(),
## county_fips_all = col_character(),
## county_name_all = col_character(),
## lat = col_double(),
## lng = col_double(),
## population = col_double(),
## density = col_double(),
## source = col_character(),
## military = col_logical(),
## incorporated = col_logical(),
## timezone = col_character(),
## ranking = col_double(),
## zips = col_character(),
## id = col_double()
## )##Question 2
###Step 2
floods = read_csv('~/github/geog-176A-labs/data/palo-flood.csv')## Parsed with column specification:
## cols(
## entityId = col_character(),
## acquisitionDate = col_datetime(format = ""),
## cloudCover = col_double(),
## processingLevel = col_character(),
## path = col_double(),
## row = col_double(),
## min_lat = col_double(),
## min_lon = col_double(),
## max_lat = col_double(),
## max_lon = col_double(),
## download_url = col_character()
## )files = lsat_scene_files(floods$download_url) %>%
filter(grepl(paste0('B',1:6,'.TIF$', collapse = '|'), file)) %>%
arrange(file) %>%
pull(file)
###Step 3
st = sapply(files, lsat_image) ## File in cache## File in cache
## File in cache
## File in cache
## File in cache
## File in cacheb = stack(st) %>% setNames(paste0('band', 1:6))
area(b)## Warning in .local(x, ...): This function is only useful for Raster* objects with
## a longitude/latitude coordinates## class : RasterLayer
## dimensions : 7811, 7681, 59996291 (nrow, ncol, ncell)
## resolution : 30, 30 (x, y)
## extent : 518085, 748515, 4506585, 4740915 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=15 +datum=WGS84 +units=m +no_defs
## source : memory
## names : layer
## values : 900, 900 (min, max)####Step 3 - What are the dimensions of your stacked image? What is the CRS? What is the cell resolution?
#### The dimensions are 7811 rows, 7861 columns, and 59996291 cells with 6 layers. The assigned CRS was WGS84 datum, UM. The cell resolution is 30 x 30.
####Step 4
cropped = palo_and_area %>%
st_as_sf() %>%
st_transform(crs(b))
bCrop = crop(b, cropped)
area(bCrop)## Warning in .local(x, ...): This function is only useful for Raster* objects with
## a longitude/latitude coordinates## class : RasterLayer
## dimensions : 340, 346, 117640 (nrow, ncol, ncell)
## resolution : 30, 30 (x, y)
## extent : 594075, 604455, 4652535, 4662735 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=15 +datum=WGS84 +units=m +no_defs
## source : memory
## names : layer
## values : 900, 900 (min, max)####The dimensions are 340 rows, 346 columns, and 117640 cells with 6 layers. The assigned CRS was WGS84, UM. The cell resolution is 30 x 30. Question 3
### Step 1
bands <- bCrop %>% setNames(c('coastal', 'blue', 'green', 'red', 'nir', 'swir'))
#### Natural Color
plotRGB(bands, r = 4, g = 3, b = 2)
#### NIR-R-G (Infared)
plotRGB(bands, r = 5, g = 4, b = 2)
#### NIR-SWIR1-R (false color)
plotRGB(bands, r = 5, g = 6, b = 4)
#### False color for vegetation and water
plotRGB(bands, r = 5, g = 7, b = 1)## Warning in .local(x, ...): layer was changed to 6
### Step 2
### Natural Color
plotRGB(bands, r = 4, g = 3, b = 2, stretch = 'lin')
#### NIR-R-G (Infared)
plotRGB(bands, r = 5, g = 4, b = 2, stretch = 'hist')
#### NIR-SWIR1-R (false color)
plotRGB(bands, r = 5, g = 6, b = 4, stretch = 'lin')
#### False color for vegetation and water
plotRGB(bands, r = 5, g = 7, b = 1, stretch = 'hist')## Warning in .local(x, ...): layer was changed to 6
#### The stretch function enhances an image by adjusting brightness, contrast, etc. For example, "lin" and "hist" are used to refer to linear or histogram stretches respectively. Different stretches may emphasize different aspects of the same plot, which will hopefully allow for ease of analysis. For example, when comparing lin and hist stretches for my last plot, false color for vegetation and water, there is a significant difference in the two. As I chose hist, the smaller waterways that split off become significantly easier to view, as the red and blue allows for a larger contrast against the green. Question 4
### Step 1
ndvi = (bands$nir - bands$red)/(bands$nir + bands$red)
ndwi = (bands$green - bands$nir)/(bands$green + bands$nir)
mndwi = (bands$green - bands$swir)/(bands$green + bands$swir)
wri = (bands$green + bands$red)/(bands$nir + bands$swir)
swi = 1/(sqrt(bands$blue-bands$swir))## Warning in sqrt(getValues(x)): NaNs producedstacks <- stack(ndvi, ndwi, mndwi, wri, swi) %>%
setNames(c('NDVI', 'NDWI', 'MNDWI', 'WRI', 'SWI'))
plot(stacks, col = colorRampPalette(c("blue", "white", "red"))(256))
#### The basic shape of the waterway is visible in all. The basic outlines of the city is visible to some degree within the first 4, but not in SWI. NDVI and NDWI have a relatively equal split of the three colors, with the water represented mostly by one color and the land areas a split of different colors, whereas MNDWI, WRI, and SWI are mostly just two colors, with the water being one color and the land being different shades, or all one shade (in the case of SWI).
### Step 2
ndviThreshold =
function(x) {ifelse(x<=0,1,0)}
ndwiThreshold =
function(x) {ifelse(x>=0,1,0)}
mndwiThreshold =
function(x) {ifelse(x>=0,1,0)}
wriThreshold =
function(x) {ifelse(x>=1,1,0)}
swiThreshold =
function(x) {ifelse(x<=5,1,0)}
ndviFlood =
calc(ndvi, ndviThreshold)
ndwiFlood =
calc(ndwi, ndwiThreshold)
mndwiFlood =
calc(mndwi, mndwiThreshold)
wriFlood =
calc(wri, wriThreshold)
swiFlood =
calc(swi, swiThreshold)
stackFlood =
stack(c(ndviFlood, ndwiFlood, mndwiFlood, wriFlood, swiFlood)) %>%
setNames(c('NDVI', 'NDWI', 'MNDWI', 'WRI', 'SWI'))
plot(stackFlood, col =
colorRampPalette(c("white","blue"))(256))
sum(is.na(values(ndviFlood)))## [1] 0sum(is.na(values(ndwiFlood)))## [1] 0sum(is.na(values(mndwiFlood)))## [1] 0sum(is.na(values(wriFlood)))## [1] 0sum(is.na(values(swiFlood)))## [1] 102433values(swiFlood) <- ifelse(is.na(values(swiFlood)), 0, values(swiFlood))Question 5
### Step 1
set.seed(20200900)
### Step 2
values <- getValues(bands)
dim(values)## [1] 117640 6#### This shows the values are expressed as a matrix of size 117640 by 6. The dimensions of bands were 340 by 346 with 6 layers, 340*246 = 117640, which means the values are stored as individual points as one dimension and the layers as the other.
values = na.omit(values)
rastClust <- kmeans(values, 12, iter.max = 100)
kmeans_raster = stacks$NDVI
values(kmeans_raster) = rastClust$cluster
plot(kmeans_raster)
### Step 3
values2 = values(ndwiFlood)
table(values2, values(kmeans_raster)) %>%
which.max()## [1] 21idx = which.max(table1[2,])## Warning in which.max(table1[2, ]): NAs introduced by coercionthresh = function(x) {ifelse(x == idx, 0, 1)}
flood = calc(kmeans_raster, thresh)
stackFlood <- addLayer(stackFlood, flood)
names(stackFlood)[6] = "K Means"
plot(stackFlood, col =
colorRampPalette(c("white","blue"))(256))
Question 6
### Step 1
kabletable = cellStats(stackFlood, sum)
knitr::kable(kabletable, caption = "Number of Flooded Cells per Image", col.names = c("Number of Cells"))| Number of Cells | |
|---|---|
| NDVI | 6666 |
| NDWI | 7212 |
| MNDWI | 11939 |
| WRI | 8469 |
| SWI | 15201 |
| K.Means | 109930 |
#### As I have answered before, the resolution is 30x30 and we know the units are meters because of the CRS, so the area is 900 meters. So a single cell has an area of 900 meters squared, so to get the total area, we need to multiply the number of cells by 900.
kabletable2 = kabletable * 900
knitr::kable(kabletable2, caption = "Area of Flooded Cells (meters squared)", col.names = c("Flooded Area"))| Flooded Area | |
|---|---|
| NDVI | 5999400 |
| NDWI | 6490800 |
| MNDWI | 10745100 |
| WRI | 7622100 |
| SWI | 13680900 |
| K.Means | 98937000 |
### Step 2
cR <-
calc(stackFlood, function(x){sum(x)}) %>%
setNames('Flood Uncertainty')
cRp = plot(cR, col = blues9, main = 'Flood Uncertainty') 
values(stackFlood) <-
ifelse(values(stackFlood)==0, NA, values(stackFlood))
### Step 3
mapview(cR)## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"): Discarded
## ellps WGS 84 in CRS definition: +proj=merc +a=6378137 +b=6378137 +lat_ts=0
## +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null +wktext +no_defs## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"): Discarded
## datum WGS_1984 in CRS definition## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"): Discarded
## ellps WGS 84 in CRS definition: +proj=merc +a=6378137 +b=6378137 +lat_ts=0
## +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null +wktext +no_defs## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"): Discarded
## datum WGS_1984 in CRS definition#### Some of the cells are not an even number, as they may